home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / raytrace / radiance / nextrad.lha / NeXtRad / vector-matrix.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-02-22  |  8.9 KB  |  307 lines
  1. /* vector-matrix.c  */
  2. /* written by : Jason R. Wilson   2/21/93 */
  3.  
  4. /* this package provides basic vector and matrix manipulations */
  5.  
  6. #include <math.h>
  7. #include <stdio.h>
  8. #include "datastruct.h"
  9.  
  10. /*****************************************************************************/
  11. /* Function:  Length                                                         */
  12. /* Purpose :  To return the length of a 3-D vector.                          */
  13. /*****************************************************************************/
  14. double Length(Vector vec) {
  15.  
  16.    return sqrt(vec.dx*vec.dx + vec.dy*vec.dy + vec.dz*vec.dz);
  17. }
  18.  
  19.  
  20. /*****************************************************************************/
  21. /* Function:  Normalize                                                      */
  22. /* Purpose :  To return the normalized version of a vector.                  */
  23. /*****************************************************************************/
  24. Vector Normalize(Vector vec) {
  25.  
  26.    Vector   temp;
  27.    double    len;
  28.  
  29.    len = Length(vec);
  30.    if (len == 0.0) 
  31.       printf("zero vector in Normalize");
  32.    else {
  33.       temp.dx = vec.dx/len;
  34.       temp.dy = vec.dy/len;
  35.       temp.dz = vec.dz/len;
  36.    }
  37.    return temp;
  38. }
  39.  
  40.  
  41. /*****************************************************************************/
  42. /* Function:  Scale                                                          */
  43. /* Purpose :  To scalar multiply the components of a 3-D vector by the       */
  44. /*            doubleing point number, scale.                                  */
  45. /*****************************************************************************/
  46. Vector Scale(Vector vec,double scale) {
  47.  
  48.    Vector  temp;
  49.  
  50.    temp.dx = scale*vec.dx;
  51.    temp.dy = scale*vec.dy;
  52.    temp.dz = scale*vec.dz;
  53.    return temp;
  54. }
  55.  
  56. /*****************************************************************************/
  57. /* Function:  Sub                                                            */
  58. /* Purpose :  To return the difference of vec1 and vec2 ==> vec1 - vec2.     */
  59. /*****************************************************************************/
  60. Vector Sub(Vector vec1, Vector vec2) {
  61.  
  62.    Vector temp;
  63.  
  64.    temp.dx = vec1.dx - vec2.dx;
  65.    temp.dy = vec1.dy - vec2.dy;
  66.    temp.dz = vec1.dz - vec2.dz;
  67.    return temp;
  68. }
  69.  
  70.  
  71. /*****************************************************************************/
  72. /* Function:  Add                                                            */
  73. /* Purpose :  To return the sum of vec1 and vec2 ==> vec1 + vec2.            */
  74. /*****************************************************************************/
  75. Vector Add(Vector vec1, Vector vec2) {
  76.  
  77.    Vector temp;
  78.  
  79.    temp.dx = vec1.dx + vec2.dx;
  80.    temp.dy = vec1.dy + vec2.dy;
  81.    temp.dz = vec1.dz + vec2.dz;
  82.    return temp;
  83. }
  84.  
  85.  
  86. /*****************************************************************************/
  87. /* Function:  Dotprod                                                        */
  88. /* Purpose :  To return the dot product of vectors, vec1 and vec2.           */
  89. /*****************************************************************************/
  90. double Dotprod(Vector vec1, Vector vec2) {
  91.  
  92.    return (vec1.dx*vec2.dx + vec1.dy*vec2.dy + vec1.dz*vec2.dz);
  93. }
  94.  
  95.  
  96. /*****************************************************************************/
  97. /* Function:  Cross                                                          */
  98. /* Purpose :  To return the cross product of vectors, vec1 and vec2.         */
  99. /*****************************************************************************/
  100. Vector Cross(Vector vec1, Vector vec2) {
  101.  
  102.    Vector temp;
  103.  
  104.    temp.dx = vec1.dy*vec2.dz - vec1.dz*vec2.dy;
  105.    temp.dy = vec1.dz*vec2.dx - vec1.dx*vec2.dz;
  106.    temp.dz = vec1.dx*vec2.dy - vec1.dy*vec2.dx;
  107.    return temp;
  108. }
  109.  
  110. /***********************************************************************/
  111. /* Function: Transpose                                                 */
  112. /* this routine puts the transpose of matrix a into matrix b           */
  113. /***********************************************************************/
  114. void Transpose (Matrix *a,Matrix *b) {
  115.  
  116.   b->mat[0] = a->mat[0];
  117.   b->mat[1] = a->mat[4];
  118.   b->mat[2] = a->mat[8];
  119.   b->mat[3] = a->mat[12];
  120.   b->mat[4] = a->mat[1];
  121.   b->mat[5] = a->mat[5];
  122.   b->mat[6] = a->mat[9];
  123.   b->mat[7] = a->mat[13];
  124.   b->mat[8] = a->mat[2];
  125.   b->mat[9] = a->mat[6];
  126.   b->mat[10] = a->mat[10];
  127.   b->mat[11] = a->mat[14];
  128.   b->mat[12] = a->mat[3];
  129.   b->mat[13] = a->mat[7];
  130.   b->mat[14] = a->mat[11];
  131.   b->mat[15] = a->mat[15];
  132.  
  133. }
  134. /*****************************************************************************/
  135. /* Function:  pointMultiply                                                  */
  136. /* Purpose :  To vector multiply the n by n matrix, b, by the n-element      */
  137. /*            vector, a, to produce the n-element vector c.                  */
  138. /*****************************************************************************/
  139. void pointMultiply(Point a, double *b, double *c) {
  140.  
  141.    c[0] = a.x * b[0] + a.y * b[4] + a.z * b[8] + 1.0 * b[12];
  142.    c[1] = a.x * b[1] + a.y * b[5] + a.z * b[9] + 1.0 * b[13];
  143.    c[2] = a.x * b[2] + a.y * b[6] + a.z * b[10] + 1.0 * b[14];
  144.    c[3] = a.x * b[3] + a.y * b[7] + a.z * b[11] + 1.0 * b[15];
  145. }
  146.  
  147. /***********************************************************************/
  148. void HomogMultiply (double *a,double *b)
  149. /* multiplies a homog. point by a 4X4 matrix */
  150. {
  151.   double temp[4];
  152.   temp[0] = a[0] * b[0] + a[1] * b[4] + a[2] * b[8] + a[3] * b[12];
  153.   temp[1] = a[0] * b[1] + a[1] * b[5] + a[2] * b[9] + a[3] * b[13];
  154.   temp[2] = a[0] * b[2] + a[1] * b[6] + a[2] * b[10] + a[3] * b[14];
  155.   temp[3] = a[0] * b[3] + a[1] * b[7] + a[2] * b[11] + a[3] * b[15];
  156.   a[0] = temp[0];
  157.   a[1] = temp[1];
  158.   a[2] = temp[2];
  159.   a[3] = temp[3];
  160. }
  161.   
  162. /***********************************************************************/
  163.   
  164. void matrixMultiply (int n, double *a, double *b, double *c)
  165. /* Multiplies the nxn matrix a by the nxn matrix b to produce the nxn matrix
  166.    c. */
  167. {
  168.    int i, j, k, posa, posb, posc;
  169.    double sum;
  170.    
  171.    for (i=0; i<n; i++)
  172.       for (j=0; j<n; j++)
  173.       {
  174.          posc = i*n + j;
  175.          sum = 0.0;
  176.      for (k=0; k<n; k++)
  177.      {
  178.         posa = i*n + k;
  179.         posb = k*n + j;
  180.         sum += a[posa] * b[posb];
  181.      }
  182.      c[posc] = sum;
  183.       }
  184. }
  185.  
  186.  
  187. void createIdentity (int n, double *a)
  188. /* Creats an nxn identity matrix in a. */
  189. {
  190.    int i, j;
  191.    
  192.    for (i=0; i<n; i++)
  193.       for (j=0; j<n; j++)
  194.          a[i*n + j] = 0.0;
  195.    for (i=0; i<n; i++)
  196.       a[i*n + i] = 1.0;
  197. }
  198.  
  199. void matrixCopy (double *a, double *b)
  200. /* copies matrix a into matrix b (assumes 4X4) */
  201. {
  202.    int loop;
  203.  
  204.    for (loop = 0;loop < 16;loop++)
  205.       b[loop] = a[loop];
  206. }
  207.  
  208. void ProduceTransform (Point VRP,Vector u, Vector v, Vector n,
  209.                double Eye,Matrix *Composition)
  210.  
  211. /* build the transform matrix that transforms coordinate systems and
  212.    does perspective projects */
  213.  
  214. {
  215.    Matrix ViewTransform, Perspective;
  216.    Vector rprime, r;
  217.  
  218.    r.dx = VRP.x;
  219.    r.dy = VRP.y;
  220.    r.dz = VRP.z;
  221.  
  222.    rprime.dx = (-1)*Dotprod(r,u);
  223.    rprime.dy = (-1)*Dotprod(r,v);
  224.    rprime.dz = (-1)*Dotprod(r,n);
  225.  
  226.    createIdentity (4,ViewTransform.mat);
  227.  
  228.    ViewTransform.mat[0] = u.dx;
  229.    ViewTransform.mat[1] = v.dx;
  230.    ViewTransform.mat[2] = n.dx;
  231.    ViewTransform.mat[4] = u.dy;
  232.    ViewTransform.mat[5] = v.dy;
  233.    ViewTransform.mat[6] = n.dy;
  234.    ViewTransform.mat[8] = u.dz;
  235.    ViewTransform.mat[9] = v.dz;
  236.    ViewTransform.mat[10] = n.dz;
  237.    ViewTransform.mat[12] = rprime.dx;
  238.    ViewTransform.mat[13] = rprime.dy;
  239.    ViewTransform.mat[14] = rprime.dz;
  240.  
  241.    /* Build the perspective transform matrix */
  242.  
  243.    createIdentity (4,Perspective.mat);
  244.    Perspective.mat[11] = (-1)/Eye;
  245.    
  246.    /* Compose the 2 matrices */
  247.  
  248.    matrixMultiply (4,ViewTransform.mat,Perspective.mat,Composition->mat);
  249.  
  250. }
  251.  
  252.  
  253. /*----------------------------------------------------------------------*/
  254.  
  255. void TransformView (ObjectCell *ObjectHead,Matrix Composition)
  256. /* This routine transforms every vertex into view coordinates, */
  257. /* and does perspective transformation */
  258.  
  259. {
  260.    VertexCell *CurrentVert;
  261.    ObjectCell *CurrentObject;
  262.  
  263.    CurrentObject = ObjectHead;
  264.  
  265.    while (!(CurrentObject == NULL))
  266.    {
  267.       CurrentVert = CurrentObject->VertexHead;
  268.       
  269.       while (!(CurrentVert == NULL))
  270.       {
  271.      pointMultiply (CurrentVert->WorldPosition,Composition.mat,
  272.             CurrentVert->ViewPosition.hom);
  273.      CurrentVert = CurrentVert->Next;
  274.       } 
  275.       CurrentObject = CurrentObject->Next;
  276.    }
  277. }
  278.  
  279.  
  280. /*----------------------------------------------------------------------*/
  281.  
  282.  
  283. void TransformObject (ObjectCell *Object)
  284. /* This routine transforms each vertex of the object from object coords. 
  285.    to world coordinates */
  286.  
  287. {
  288.    VertexCell *CurrentVert;
  289.    double tempHom[4];
  290.  
  291.    CurrentVert = Object->VertexHead;
  292.    
  293.    while (!(CurrentVert == NULL))
  294.    {
  295.       pointMultiply (CurrentVert->WorldPosition,Object->Transform.mat,tempHom);
  296.       CurrentVert->WorldPosition.x = tempHom[0];
  297.       CurrentVert->WorldPosition.y = tempHom[1];
  298.       CurrentVert->WorldPosition.z = tempHom[2];
  299.       CurrentVert = CurrentVert->Next;
  300.    } 
  301. }
  302.  
  303.  
  304. /***********************************************************************/
  305.  
  306.  
  307.